dfp <- df.eel %>%filter(fi_lfs_code =="S", fi_year >1997,# norway has three silver eels in total, so I remove this!ecoreg =="norwegian") %>%rename(length = lengthmm) %>%mutate(sex =as.integer(if_else(sex =="f", 0, 1)))
filter: removed 724,965 rows (83%), 147,494 rows remaining
rename: renamed one variable (length)
mutate: converted 'sex' from character to integer (0 new NA)
mutate: new variable 'mb' (integer) with 154 unique values and 0% NA
new variable 'er' (integer) with 9 unique values and 0% NA
new variable 'year' (integer) with 28 unique values and 0% NA
new variable 'temp.sc' (double) with 78 unique values and 0% NA
new variable 'age.sc' (double) with one unique value and 100% NA
filter: removed 131,289 rows (89%), 16,205 rows remaining
mutate: converted 'sex' from integer to character (0 new NA)
Show code
# possible temperature effects on silvering lengthdata.silv %>%mutate(sex =if_else(sex ==0, "f", "m")) %>%ggplot(aes(tmp_dc_uyr/10, length)) +geom_point(show.legend =FALSE, alpha =0.5) +geom_smooth(method ="lm", se =TRUE) +facet_wrap(~factor(sex)) +theme_minimal() +labs(x ="Temp. [C]")
mutate: converted 'sex' from integer to character (0 new NA)
`geom_smooth()` using formula = 'y ~ x'
Source models / Load samples
Show code
# Load samples# readRDS(paste0(home,"/data/silv.samples_1214.RData"))$WAIC# readRDS(paste0(home,"/data/silv.samples_1216.RData"))$WAICsilv.samples <-readRDS(paste0(home,"/models_eel/samples/silv.samples_2026-01-08.RData"))$samples # v9# # Or source model code and mcmc (load libs and data)# t <- Sys.time()# source(file = paste0(home,"/Eel models/2-2_model_eel_silver_v5.R"))# Sys.time() - t # approx xxx hours with xx million hmc samples
Traces, ess, and ‘potential scale reduction factor’
node.sub <-grep("^R\\[", colnames(silv.samples[[1]]), value =TRUE)silv.samples[, node.sub[sample(1:length(node.sub), 36)], drop =FALSE] %>%mcmc_trace()
Show code
pl <- silv.samples[, node.sub[], drop =FALSE] %>%gelman.diag(multivariate =FALSE) # Error in chol.default(W) : the leading minor of order 1 is not positivepl <- silv.samples[, node.sub[], drop =FALSE] %>%effectiveSize()print("effective sample size")
mutate (grouped): new variable 'male' (double) with 377,496 unique values and 0% NA
new variable 'year' (double) with 28 unique values and 0% NA
rename: renamed one variable (female)
pivot_longer: reorganized (female, male) into (sex, length) [was 378000x9, now 756000x9]
mutate (grouped): changed 756,000 values (100%) of 'length' (0 new NAs)
distinct: removed 147,485 rows (>99%), 9 rows remaining
Joining with `by = join_by(er)`
left_join: added one column (ecoreg)
> rows only in x 0
> rows only in data.silv %>% distinct(.. ( 0)
> matched rows 756,000
> =========
> rows total 756,000
ungroup: no grouping variables remain
mutate: new variable 'm_length' (double) with 56 unique values and 0% NA
Show code
dfp %>%filter(sex =="female") %>%ungroup() %>%ggplot() +geom_boxplot(aes(y = length, x = year, group = year), color ="#F8766D", outliers =FALSE) +geom_line(aes(y = m_length, x = year), color ="#F8766D", alpha =0.4) +facet_wrap(~ecoreg, scales ="free_y") +theme_minimal() +theme(axis.text.x =element_text(angle =90, hjust =1, size =rel(.75))) +labs(x ="year")
ggsave(paste0(home, "/report/silv_lm.png"), width =17, height =14, units ="cm")silv.samples %>%spread_draws(alpha[er, yr], bs[er],sep =",") %>%median_qi() %>%mutate(male = alpha + bs,year = yr +1997) %>%rename(female = alpha) %>%pivot_longer(cols =c("female","male"), names_to ="sex", values_to ="length") %>%left_join(data.silv %>%distinct(er,ecoreg)) %>%mutate(length_mc =exp(length) -mean(exp(length)), .by =c(sex)) %>%ggplot() +geom_line(aes(y = length_mc, x = year, color = ecoreg)) +facet_wrap(~sex, scales ="free_y") +theme_minimal() +labs(x ="year", caption ="centered by mean by sex")
mutate: new variable 'male' (double) with 252 unique values and 0% NA
new variable 'year' (double) with 28 unique values and 0% NA
rename: renamed one variable (female)
pivot_longer: reorganized (female, male) into (sex, length) [was 252x13, now 504x13]
distinct: removed 147,485 rows (>99%), 9 rows remaining
Joining with `by = join_by(er)`
left_join: added one column (ecoreg)
> rows only in x 0
> rows only in data.silv %>% distinct(.. ( 0)
> matched rows 504
> =====
> rows total 504
mutate: new variable 'length_mc' (double) with 504 unique values and 0% NA
mutate: new variable 'r' (integer) with 9 unique values and 0% NA
pivot_longer: reorganized (V1, V2, V3, V4, V5, …) into (c, value) [was 9x10, now 81x3]
mutate: converted 'c' from character to integer (0 new NA)
left_join: added 2 columns (ecoreg, m.lat)
> rows only in x 0
> rows only in var.lats ( 0)
> matched rows 81
> ====
> rows total 81
left_join: added 4 columns (ecoreg.x, m.lat.x, ecoreg.y, m.lat.y)
> rows only in x 0
> rows only in var.lats ( 0)
> matched rows 81
> ====
> rows total 81
mutate: converted 'r' from integer to factor (0 new NA)
converted 'c' from integer to factor (0 new NA)
